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We address the problem of reconstructing a multi-band signal from its sub-Nyquist point-wise samples. To 
date, all reconstruction methods proposed for this class of signals assumed knowledge of the band locations. In this 
paper, we develop a non-linear blind perfect reconstruction scheme for multi-band signals which does not require 
the band locations. Our approach assumes an existing blind multi-coset sampling method. The sparse structure 
of multi-band signals in the continuous frequency domain is used to replace the continuous reconstruction with 
a single finite dimensional problem without the need for discretization. The resulting problem can be formulated 
within the framework of compressed sensing, and thus can be solved efficiently using known tractable algorithms 
from this emerging area. We also develop a theoretical lower bound on the average sampling rate required for 
blind signal reconstruction, which is twice the minimal rate of known-spectrum recovery. Our method ensures 
^sO , perfect reconstruction for a wide class of signals sampled at the minimal rate. Numerical experiments are presented 

demonstrating blind sampling and reconstruction with minimal sampling rate. 
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Index Terms 

Kruskal-rank, Landau-Nyquist rate, multiband, multiple measurement vectors (MMV), nonuniform periodic 
sampling, orthogonal matching pursuit (OMP), signal representation, sparsity. 

I. Introduction 

The well known Whittaker, Kotelriikov, and Shannon (WKS) theorem links analog signals with a discrete 
representation, allowing the transfer of the signal processing to a digital framework. The theorem states that a 
real-valued signal bandlimited to B Hertz can be perfectly reconstructed from its uniform samples if the sampling 
rate is at least 2B samples per second. This minimal rate is called the Nyquist rate of the signal. 

Multi-band signals are bandlimited signals that posses an additional structure in the frequency domain. The 
spectral support of a multi-band signal is restricted to several continuous intervals. Each of these intervals is called 
a band and it is assumed that no information resides outside the bands. The design of sampling and reconstruction 
systems for these signals involves three major considerations. One is the sampling rate. The other is the set of 
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multi-band signals that the system can perfectly reconstruct. The last one is blindness, namely a design that does 
not assume knowledge of the band locations. Blindness is a desirable property as signals with different band 
locations are processed in the same way. Landau [1] developed a minimal sampling rate for an arbitrary sampling 
method that allows perfect reconstruction. For multi-band signals, the Landau rate is the sum of the band widths, 
which is below the corresponding Nyquist rate. 

Uniform sampling of a real bandpass signal with a total width of 2B Hertz on both sides of the spectrum was 
studied in [2], It was shown that only special cases of bandpass signals can be perfectly reconstructed from their 
uniform samples at the minimal rate of 2B samples/sec. Kohlenberg [3] suggested periodic non-uniform sampling 
with an average sampling rate of 2B. He also provided a reconstruction scheme that recovers any bandpass 
signal exactly. Lin and Vaidyanathan [4] extended his work to multi-band signals. Their method ensures perfect 
reconstruction from periodic non uniform sampling with an average sampling rate equal to the Landau rate. Both 
of these works lack the blindness property as the information about the band locations is used in the design of 
both the sampling and the reconstruction stages. 

Herley and Wong [5] and Venkataramani and Bresler [8] suggested a blind multi-coset sampling strategy that 
is called universal in [8]. The authors of [8] also developed a detailed reconstruction scheme for this sampling 
strategy, which is not blind as its design requires information about the spectral support of the signal. Blind 
multi-coset sampling renders the reconstruction applicable to a wide set of multi-band signals but not to all of 
them. 

Although spectrum-blind reconstruction was mentioned in two conference papers in 1996 [6], [7], a full spectrum- 
blind reconstruction scheme was not developed in these papers. It appears that spectrum-blind reconstruction has 
not been handled since then. 

We begin by developing a lower bound on the minimal sampling rate required for blind perfect reconstruction 
with arbitrary sampling and reconstruction. As we show the lower bound is twice the Landau rate and no more 
than the Nyquist rate. This result is based on recent work of Lue and Do [20] on sampling signals from a union 
of subspaces. 

The heart of this paper is the development of a spectrum-blind reconstruction (SBR) scheme for multi-band 
signals. We assume a blind multi-coset sampling satisfying the minimal rate requirement. Theoretical tools are 
developed in order to transform the continuous nature of the reconstruction problem into a finite dimensional 
problem without any discretization. We then prove that the solution can be obtained by finding the unique sparsest 
solution matrix from Multiple-Measurement- Vectors (MMV). This set of operations is grouped under a block we 
name Continuous to Finite (CTF). This block is the cornerstone of two SBR algorithms we develop to reconstruct 
the signal. One is entitled SBR4 and enables perfect reconstruction using only one instance of the CTF block 
but requires twice the minimal sampling rate. The other is referred to as SBR2 and allows for sampling at the 
minimal rate, but involves a bi-section process and several uses of the CTF block. Other differences between the 
algorithms are also discussed. Both SBR4 and SBR2 can easily be implemented in DSP processors or in software 
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environments. 

Our proposed reconstruction approach is applicable to a broad class of multi-band signals. This class is the blind 
version of the set of signals considered in [8]. In particular, we characterize a subset M. of this class by the maximal 
number of bands and the width of the widest band. We then show how to choose the parameters of the multi-coset 
stage so that perfect reconstruction is possible for every signal in M.. This parameter selection is also valid for 
known-spectrum reconstruction with half the sampling rate. The set M. represents a natural characterization of 
multi-band signals based on their intrinsic parameters which are usually known in advance. We prove that the 
SBR4 algorithm ensures perfect reconstruction for all signals in A4. The SBR2 approach works for almost all 
signals in M. but may fail in some very special cases (which typically will not occur). As our strategy is applicable 
also for signals that do not lie in M., we present a nice feature of a success recovery indication. Thus, if a signal 
cannot be recovered this indication prevents further processing of invalid data. 

The CTF block requires finding a sparsest solution matrix which is an NP-hard problem [12]. Several sub-optimal 
efficient methods have been developed for this problem in the compressed sensing (CS) literature [15], [16]. In our 
algorithms, any of these techniques can be used. Numerical experiments on random constructions of multi-band 
signals show that both SBR4 and SBR2 maintain a satisfactory exact recovery rate when the average sampling 
rate approaches their theoretical minimum rate requirement and sub-optimal implementations of the CTF block 
are used. Moreover, the average runtime is shown to be fast enough for practical usage. 

Our work differs from other main stream CS papers in two aspects. The first is that we aim to recover a 
continuous signal, while the classical problem addressed in the CS literature is the recovery of discrete and finite 
vectors. An adaptation of CS results to continuous signals was also considered in a set of conferences papers (see 
[21], [22] and the references therein). However, these papers did not address the case of multi-band signals. In 
[22] an underlying discrete model was assumed so that the signal is a linear combination of a finite number of 
known functions. Here, there is no discrete model as the signals are treated in a continuous framework without any 
discretization. The second aspect is that we assume a deterministic sampling stage and our theorems and results 
do not involve any probability model. In contrast, the common approach in compressive sensing assumes random 
sampling operators and typical results are valid with some probability less than 1 [13], [19], [21], [22]. 

The paper is organized as follows. In Section II we formulate our reconstruction problem. The minimal density 
theorem for blind reconstruction is stated and proved in Section III. A brief overview of multi-coset sampling is 
presented in Section IV. We develop our main theoretical results on spectrum-blind reconstruction and present the 
CTF block in Section V. Based on these results, in Section VI, we design and compare the SBR4 and the SBR2 
algorithms. Numerical experiments are described in Section VII. 
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II. Preliminaries and Problem formulation 

A. Notation 

Common notation, as summarized in Table I, is used throughout the paper. Exceptions to this notation are 
indicated in the text. 

TABLE I 
Notation 



x(t) 


continuous time signal with finite energy 


x(f) 


Fourier transform of x(t) (that is assumed to exist) 


a[n] 


bounded energy sequence 


z* 


conjugate of the complex number z 


V 


vector 


Vj or v(i) 


ith entry of v 


v(/) 


vector that depends on a continuous parameter / 


A 


matrix 




ikth entry of A 


A T ,A H 


transpose and the conjugate-transpose of A 


A y o 


A is an Hermitian positive semi-definite (PSD) matrix 


At 


the Moore-Penrose pseudo-inverse of A 


S 


finite or countable set 


Si 


ith element of S 


\s\ 


cardinality of a finite set S 


T 


infinite non-countable set 


A(T) 


the Lebesgue measure of T C R 



In addition, the following abbreviations are used. The t p norm of a vector v is defined as 

l|v||£ = EN P > p^°- 

i 

The default value for p is 2, so that ||v|| denotes the £2 norm of v. The standard L2 norm is used for continuous 
signals. The ith column of A is written as Aj, the ith row is (A T )j written as a column vector. 
Indicator sets for vectors and matrices are defined respectively as 

I(v) = {k I v(fc) + 0}, 1(A) = {k I (A T ) k + 0}. 

The set /(v) contains the indices of non-zero values in the vector v. The set 1(A) contains the indices of the 
non-identical zero rows of A. 

Finally, Ag is the matrix that contains the columns of A with indices belonging to the set S. The matrix Ag 
is referred to as the (columns) restriction of A to S. Formally, 

(A s ) i = (A) Si , 1<*<|5|. 
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Similarly, A is referred to as the rows restriction of A to S. 
B. Multi-band signals 

In this work our prime focus is on the set A4 of all complex-valued multi-band signals bandlimited to T = 
[0, l/T] with no more than N bands where each of the band widths is upper bounded by B. Fig. 1 depicts a 
typical spectral support for x(t) G M. 

B 
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Fig. 1. Typical spectrum support of x(t) £». 

The Nyquist rate corresponding to any x(t) € M. is l/T. The Fourier transform of a multi-band signal has 
support on a finite union of disjoint intervals in T . Each interval is called a band and is uniquely represented by 
its edges [aj,&j]. Without loss of generality it is assumed that the bands are not overlapping. 

Although our interest is mainly in signals x(t) € M., our results are applicable to a broader class of signals, 
as explained in the relevant sections. In addition, the results of the paper are easily adopted to real-valued signals 
supported on [— 1/2T, +1/2T]. The required modifications are explained in Appendix A and are based on the 
equations derived in Section IV-A. 

C. Problem formulation 

We wish to perfectly reconstruct x(t) G M from its point-wise samples under two constraints. One is blindness, 
so that the information about the band locations is not used while acquiring the samples and neither can it be 
used in the reconstruction process. The other is that the sampling rate required to guarantee perfect reconstruction 
should be minimal. 

This problem is solved if either of its constraints is removed. Without the rate constraint, the WKS theorem 
allows perfect blind-reconstruction for every signal x(t) bandlimited to T from its uniform samples at the Nyquist 
rate x(t = n/T). Alternatively, if the exact number of bands and their locations are known, then the method of 
[4] allows perfect reconstruction for every multi-band signal at the minimal sampling rate provided by Landau's 
theorem [1]. 

In this paper, we first develop the minimal sampling rate required for blind reconstruction. We then use a multi- 
coset sampling strategy to acquire the samples at an average sampling rate satisfying the minimal requirement. The 
design of this sampling method does not require knowledge of the band locations. We provide a spectrum-blind 
reconstruction scheme for this sampling strategy in the form of two different algorithms, named SBR4 and SBR2. It 
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is shown that if the sampling rate is twice the minimal rate then algorithm SBR4 guarantees perfect reconstruction 
for every x(t) € M.. The SBR2 algorithm requires the minimal sampling rate and guarantees perfect reconstruction 
for most signals in M.. However, some special signals from M., discussed in Section VI-B, cannot be perfectly 
reconstructed by this approach. Excluding these special cases, our proposed method satisfies both constraints of 
the problem formulation. 

III. Minimal sampling rate 

We begin by quoting Landau's theorem for the minimal sampling rate of an arbitrary sampling method that 
allows known-spectrum perfect reconstruction. It is then proved that blind perfect-reconstruction requires a minimal 
sampling rate that is twice the Landau rate. 

A. Known spectrum support 

Consider the space of bandlimited functions restricted to a known support T C T: 

B T = {x{t) GL 2 (R)| sup P X(/) CT}. (1) 

A classical sampling scheme takes the values of x(t) on a known countable set of locations R = {r n }'^L- o- The 
set R is called a sampling set for Br if x(t) can be perfectly reconstructed in a stable way from the sequence of 
samples xr[ti] = x(t = r n ). The stability constraint requires the existence of constants a > and (3 < oo such 
that: 

a\\x - y\\ 2 < \\x R - y R \\ 2 < (3\\x - y\\ 2 , Vx,y€#T- (2) 
Landau [1] proved that if R is a sampling set for Br then it must have a density D~(R) > A(T), where 

D~(R)= liminfl^M±^ (3) 



is the lower Beurling density, and A(T) is the Lebesgue measure of T. The numerator in (3) counts the number 
of points from R in every interval of width r of the real axis 1 . This result is usually interpreted as a minimal 
average sampling rate requirement for Br, and A(T) is called the Landau rate. 

B. Unknown spectrum support 

Consider the set An of signals bandlimited to T with bandwidth occupation no more than < < 1, so that 

A(su PP X(/)) < Vx(t) e Nq. 



'The numerator is not necessarily finite but as the sampling set is countable the infimum takes on a finite value. 



The Nyquist rate for Afn is l/T. Note that Afn is not a subspace so that the Landau theorem is not valid here. 
Nevertheless, it is intuitive to argue that the minimal sampling rate for Afn cannot be below ft/T as this value is 
the Landau rate had the spectrum support been known. 

A blind sampling set R for Afn is a sampling set whose design does not assume knowledge of supp X(f). 
Similarly to (2) the stability of R requires the existence of a > and f3 < oo such that: 



a\\x - y\\ 2 < \\x R - y R \\ 2 < 0\\x - y\\ 2 , Vx,y G Afn- 

Theorem 1 (Minimal sampling rate): Let R be a blind sampling set for Afn- Then, 

2ft 1 



D~(R) > min 



Proof: The set Afn is of the form 



T ' T 



Mn=\J B T , 



Ter 



where 



r = {T | T C T, A(T) < JtyT}. 



(4) 

(5) 
(6) 

(7) 



Clearly, Afn is a non-countable union of subspaces. Sampling signals that lie in a union of subspaces has been 
recently treated in [20]. For every 7, 9 E T define the subspaces 



B 



7,0 — ^7 



B 1 + B e = {x + y I x e B 1 , y G 



(8) 



Since R is a sampling set for Afn, (4) holds for some constants a > 0, f3 < 00. It was proved in [20, Proposition 
2] that (4) is valid if and only if 



a\\x - y\\ 2 < \\x R - y R \\ 2 < (3\\x - y\\ 2 , Mx,y G B lfi 



(9) 



holds for every 7, 9 G T. In particular, R is a sampling set for every £> 7 # with 7, # G T. 

Observe that the space £? 7 # is of the form (1) with T = 7 U 9. Applying Landau's density theorem for each 
7, 9 G T results in 

D~(R) > A(7 U 9), V 7 ,#Gr. (10) 



Choosing 



we have that for Q, < 0.5, 



If Q > 0.5 then 7 U 9 = T and 



7 



" ft" 




1" 




, 9 = 











£>~(it:) > A( 7 U 9) = A( 7 ) + A(0) 



2ft 
I 7 ' 



D~(R) > A(7U 



1 

T' 



(11) 



(12) 
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Combining (11) and (12) completes the proof. ■ 
In [20], the authors consider minimal sampling requirements for a union of shift-invariant subspaces, with a 
particular structure of sampling functions. Specifically, they view the samples as inner products with sampling 
functions of the form {ipk(t — m)}i<k<K,m.ez> which includes multi-coset sampling. Theorem 1 extends this result 
to an arbitrary point-wise sampling operator. In particular, it is valid for non periodic sampling sets that are not 
covered by [20]. 

An immediate corollary of Theorem 1 is that if Q > 0.5 then uniform sampling at the Nyquist rate with an 
ideal low pass filter satisfies the requirements of our problem formulation. Namely, both the sampling and the 
reconstruction do not use the information about the band locations, and the sampling rate is minimal according to 
Theorem 1. As M. is contained in the space of bandlimited signals, this choice also provides perfect reconstruction 
for every x(t) € M. Therefore, in the sequel we assume that U < 0.5 so that the minimal sampling rate of 
Theorem 1 is exactly twice the Landau rate. 

It is easy to see that M. C A/q for Q, = NBT. Therefore, for known spectral support, the Landau rate is NB. 
Despite the fact that M. is a true subset of Mnbt, the proof of Theorem 1 can be adopted to show that a minimal 
density of 2NB is required so that stable perfect reconstruction is possible for signals from M.. 

We point out that both Landau's and Theorem 1 state a lower bound but do not provide a method to achieve 
the bounds. The rest of the paper is devoted to developing a reconstruction method that approaches the minimal 
sampling rate of Theorem 1 . 

IV. Universal Sampling 

This section reviews multi-coset sampling which is used in our development. We also briefly explain the 
fundamentals of known-spectrum reconstruction as derived in [8]. 

A. Multi-coset sampling 

Uniform sampling of x(t) at the Nyquist rate results in samples x(t = nT) that contain all the information 
about x{t). Multi-coset sampling is a selection of certain samples from this grid. The uniform grid is divided into 
blocks of L consecutive samples. A constant set C of length p describes the indices of p samples that are kept in 
each block while the rest are zeroed out. The set C = {cj}^ =1 is referred to as the sampling pattern where 

< ci < c 2 < ... < Cp < L - 1. (13) 

Define the ith sampling sequence for 1 < i < p as 

{x(t = nT) n = mL + for some m € Z 
otherwise. (14) 

The sampling stage is implemented by p uniform sampling sequences with period 1/(LT), where the ith sampling 
sequence is shifted by c{T from the origin. Therefore, a multi-coset system is uniquely characterized by the 



parameters L,p and the sampling pattern C. 
Direct calculations show that [8] 

L-l 



X^ T ) = ^ £ exp (j^cA X(f + { f ), (15) 



r=0 



0,^ ] 



where X c .(e^ 2n ^ T ) is the discrete-time Fourier transform (DTFT) of x Ci [n]. Thus, the goal is to choose parameters 
L,p, C such that X(f) can be recovered from (15). 

For our purposes it is convenient to express (15) in a matrix form as 

y(/) = Ax(/), V/G^b, (16) 

where y(/) is a vector of length p whose ith element is X c% {e : ' 2lT * T ), and the vector x(/) contains L unknowns 
for each / 

Xi(/) = X (f + -^V 0<i<L-l, /G^o- (17) 
The matrix A depends on the parameters L,p and the set C but not on x(t) and is defined by 

A ik = exp (j-jr-Cik) ■ (18) 

Dealing with real- valued multi-band signals requires simple modifications to (16). These adjustments are detailed 
in Appendix A. 

The Beurling lower density (i.e. the average sampling rate) of a multi-coset sampling set is 

tJ— = -pf;, (19) 

J- AVG Li 1 

which is lower than the Nyquist rate for p < L. However, an average sampling rate above the Landau rate is not 
sufficient for known-spectrum reconstruction. Additional conditions are needed as explained in the next section. 

B. Known-spectrum reconstruction and universality 

The presentation of the reconstruction is simplified using CS sparsity notation. A vector v is called i^T-sparse if 
the number of non-zero values in v is no greater than K. Using the 1$ pseudo-norm the sparsity of v is expressed 
as 1 1 v 1 1 o < K. We use the following definition of the Kruskal-rank of a matrix [14]: 

Definition 1: The Kruskal-rank of A, denoted as <r(A), is the maximal number q such that every set of q 
columns of A is linearly independent. 

Observe that for every / G the system of (16) has less equations than unknowns. Therefore, a prior on x(/) 
must be used to allow for recovery. In [8] it is assumed that the information about the band locations is available 
in the reconstruction stage. This information supplies the set 7(x(/)) for every / G Tq. Without any additional 
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prior the following condition is necessary for known-spectrum perfect reconstruction 

x(/) is p-sparse , V/ G J-q. (20) 

Using the Kruskal-rank of A a sufficient condition is formulated as 

x(/) is (r(A)-sparse , V/ € JF . (21) 

The known-spectrum reconstruction of [8] basically restricts the columns of A to I(x(/)) and inverts the resulting 
matrix in order to recover x(/). 

A sampling pattern C that yields a fully Kruskal-rank A is called universal and corresponds to <t(A) = p. 
Therefore, the set of signals that are consistent with (21) is the broadest possible if a universal sampling pattern 
is used. As we show later, choosing L < p > N and a universal pattern C makes (21) valid for every signal 
x(t) e M. 

Finding a universal pattern C, namely one that results in a fully Kruskal-rank A, is a combinatorial process. 
Several specific constructions of sampling patterns that are proved to be universal are given in [8], [10]. In particular, 
choosing L to be prime renders every pattern universal [10]. 

To summarize, choosing a universal pattern allows recovery of any x(t) satisfying (20) when the band locations 
are known in the reconstruction. We next consider blind signal recovery using universal sampling patterns. 

V. Spectrum-Blind Reconstruction 

In this section we develop the theory needed for SBR. These results are then used in the next section to construct 
two efficient algorithms for blind signal reconstruction. 

The theoretical results are devoted in the following steps: We first note that when considering blind-reconstruction, 
we cannot use the prior of [8]. In Section V-A we present a different prior that does not assume the information 
about the band locations. Using this prior we develop a sufficient condition for blind perfect reconstruction which 
is very similar to (21). Furthermore, we prove that under certain conditions on L,p,C, perfect reconstruction is 
possible for every signal in M.. We then present the basic SBR paradigm in Section V-B. The main result of the 
paper is transforming the continuous system of (16) into a finite dimensional problem without using discretization. 
In Section V-C we develop two propositions for this purpose, and present the CTF block. 

A. Conditions for blind perfect reconstruction 

Recall that for every / € Tq the system of (16) is undetermined since there are fewer equations than unknowns. 
The prior assumed in this paper is that for every / 6 Tq the vector x(/) is sparse but in contrast to [8] the location 
of the non-zero values is unknown. Clearly, in this case (20) is still necessary for blind perfect reconstruction. The 
following theorem from the CS literature is used to provide a sufficient condition. 
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Theorem 2: Suppose x is a solution of y = Ax. If ||x||o < <r(A)/2 then x is the unique sparsest solution of 
the system. 

Theorem 2 and its proof are given in [11], [15] with a slightly different notation of Spark(-A) instead of the 
Kruskal-rank of A. Note that the condition of the theorem is not necessary as there are examples that the sparsest 
solution x of y = Ax is unique while x > er(A)/2. 

Using Theorem 2, it is evident that perfect reconstruction is possible for every signal satisfying 

x(/) is -sparse , V/ G F Q . (22) 

As before, choosing a universal pattern makes the set of signals that conform with (22) the widest possible. Note 
that a factor of two distinguishes between the sufficient conditions of (21) and of (22), and results from the fact 
that here we do not know the locations of the non-zero values in x(/). 

Note that (22) provides a condition under which perfect reconstruction is possible, however, it is still unclear 
how to find the original signal. Although the problem is similar to that described in the CS literature, here finding 
the unique sparse vector must be solved for each value / in the continuous interval JF , which clearly cannot be 
implemented. 

In practice, conditions (21) and (22) are hard to verify since they require knowledge of x(t) and depend on 
the parameters of the multi-coset sampling. We therefore prefer to develop conditions on the class M. which 
characterizes multi-band signals based on their intrinsic properties: the number of bands and their widths. It is 
more likely to know the values of N and B in advance than to know if the signals to be sampled satisfy (21) or 
(22). The following theorem describes how to choose the parameters L,p and C so that the sufficient conditions 
for perfect reconstruction hold true for every x(t) G M., namely it is a unique solution of (16). The theorem 
is valid for both known and blind reconstruction with a slight difference resulting from the factor of two in the 
sufficient conditions. 

Theorem 3 (Uniqueness): Let x(t) G M. be a multi-band signal. If: 

1) The value of L is limited by 

L < JL, (23) 

2) p > N for known reconstruction or p > 2N for blind, 

3) C is a universal pattern, 

then, for every / 6 the vector x(/) is the unique solution of (16). 

Proof: If L is limited by (23) then for the zfh band % = [en, bj\ we have 

X{%) <B<^, 1 < i < N. 

Therefore, / G % implies 
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According to (17) for every / G Tq the vector x(/) takes the values of X(f) on a set of L points spaced by 
1/LT. Consequently, the number of non-zero values in x(/) is no greater than the number of the bands, namely 
x(/) is iV-sparse. 

Since C is a universal pattern, <x(A) = p. This implies that conditions (21) and (22) are satisfied. ■ 
Note that the condition on the value of p implies the minimal sampling rate requirement. To see this, substitute 
(23) into (19): 

T7^>pB. (24) 



Tayg LT 

As pointed out in the end of Section III-B, if the signals are known to lie in M. then the Landau rate is NB, 
which is implied by p > N. Theorem 1 requires an average sampling rate of 2NB, which can be guaranteed if 

p > 2N. 

B. Reconstruction paradigm 

The goal of our reconstruction scheme is to recover the signal x(t) from the set of sequences x c% [n], 1 < i < p. 
Equivalently, the aim is to reconstruct x(/) of (16) for every / G T§ from the input data y(/). 

A straight forward approach is to find the sparsest solution x(/) on a dense grid of / G T§. However, this 
discretization strategy cannot guarantee perfect reconstruction. In contrast, our approach is exact and does not 
require discretization. 

Our reconstruction paradigm is targeted at finding the diversity set which depends on x(t) and is defined as 

s= IJ J ( x (/))- < 25 ) 

The SBR algorithms we develop in Section VI are aimed at recovering the set S. With the knowledge of S perfect 
reconstruction of x(/) is possible for every / G by noting that (16) can be written as 



If the diversity set of x(t) satisfies 



then 



y(/)=A s x s (/). (26) 



\S\ < <t(A), (27) 



(A 5 ) f A s = /, (28) 
where As is of size p x Multiplying both sides of (26) by (As)T results in: 

x 5 (/)=(A 5 )V(/), V/G^b- (29) 



From (25), 



*(/) = (), Vfef ,i<£S. 



(30) 
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Thus, once S is known, and as long as (27) holds, perfect reconstruction can be obtained by (29)-(30). 

As we shall see later on (27) is implied by the condition required to transform the problem into a finite 
dimensional one. Furthermore, the following proposition shows that for x(t) € M, (27) is implied by the parameter 
selection of Theorem 3. 

Proposition 1: If L is limited by (23) then \S\ < 2N. If in addition p > 2N and C is universal then for every 
x(t) € M, the set S satisfies (27). 

Proof: The bands are continuous intervals upper bounded by B. From (17) it follows that x(/) is constructed 
by dividing T into L equal intervals of length 1/LT. Therefore if L is limited by (23) then each band can 
either be fully contained in one of these intervals or it can be split between two consecutive intervals. Since 
the number of bands is no more than N it follows that \S\ < 2N. With the additional conditions we have that 



As we described, our general strategy is to determine the diversity set S and then recover x(t) via (29)-(30). In 
the non-blind setting, S is known, and therefore if it satisfies (27) then the same equations can be used to recover 
x(t). However, note that when the band locations are known, we may use a value of p that is smaller than 2N 
since the sampling rate can be reduced. Therefore, (27) may not hold. Nonetheless, it is shown in [8], that the 
frequency axis can be divided into intervals such that this approach can be used over each frequency interval. 
Therefore, once the set S is recovered there is no essential difference between known and blind reconstruction. 

C. Formulation of a finite dimensional problem 

The set of equations of (16) consists of an infinite number of linear systems because of the continuous variable 
/. Furthermore, the expression for the diversity set S given in (25) involves a union over the same continuous 
variable. The main result of this paper is that S can be recovered exactly using only one finite dimensional problem. 
In this section we develop the underlying theoretical results that are used for this purpose. 

Consider a given T C JF . Multiplying each side of (16) by its conjugate transpose we have 



(7(A) = p > 2N > \S 



y(/)y"(/) = Ax(/)x"(/)A 



H 



V/eT. 



(31) 



Integrating both sides over the continuous variable / gives 



Q = AZ A 



H 



(32) 



with the p x p matrix 




(33) 



and the L x L matrix 




(34) 
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Define the diversity set of the interval T as 



Now, 



S r = |J J(x(/)). (35) 

/6T 



5o)ii = / \Mf)\ 2 df. 
Jfer 



This means that (Zq)h = if and only if Xj(/) = 0, V/ G T, which implies that S7- = /(Zo). 

The next proposition is used to determine whether Zq can be found by a finite dimensional problem. The 
proposition is stated for general matrices Q, A. 

Proposition 2: Suppose Q y of size p x p and A are given matrices. Let Z be any L x L matrix satisfying 

Q = AZA H , (36a) 
Z y 0, (36b) 
|J(Z)| <<r(A). (36c) 



Then, rank(Z) = rank(Q). If, in addition, 



|/(Z)| < (36d) 



then, Z is the unique solution of (36a)-(36d). 

Proof: Let Z satisfy (36a)-(36c). Define rq = rank(Q), rz = rank(Z). Since Z y it can be decomposed 
as Z = PP^ with P of size L x rz having orthogonal columns. From (36a), 

Q = (AP)(AP)^. (37) 

It can be easily be concluded that I{Z) = I(P), and thus |/(P)| < cr(A). The following lemma whose proof is 
given in Appendix B ensures that the matrix AP of size p x rz also has full column rank. 

Lemma 1: For every two matrices A,P, if |I(P)| < cr(A) then rank(P) = rank(AP). 

Since for every matrix M it is true that rank(M) = rank(MM^), (37) implies rz = tq. 

For the second part of Proposition 2 suppose that Z, Z both satisfy (36a),(36b),(36d). From the first part, 

rank(Z) = rank(Z) = tq. 
Following the earlier decompositions we write 

Z = PP H , /(Z) = J(P) (38) 



PP ff , 7(Z) = /(P). 
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In addition, 



From (36a), 



which implies that 



|/(P)l<^, |/(P)I<^- (39) 



Q = (AP)(AP) H = (AP)(AP) H , (40) 



A(P - PR) = 0, (41) 



for some unitary matrix R. It is easy to see that (39) results in |/(PR)| < <r(A)/2. Therefore the matrix P — PR 
has at most c(A) non-identical zero rows. Applying Lemma 1 to (41) results in P = PR. Substituting this into 
(38) we have that Z = Z. ■ 
The following proposition shows how to construct the matrix Z by finding the sparsest solution of a linear 
system. 

Proposition 3: Consider the setting of Proposition 2 and assume Z satisfies (36d). Let r = rank(Q) and define 
a matrix V of size p x r using the decomposition Q = VV^, such that V has r orthogonal columns. Then the 
linear system 

V = AU (42) 

has a unique sparsest solution matrix Uo- Namely, V = AUo and |/(Uo)| is minimal. Moreover, Z = UoU^. 

Proof: Substitute the decomposition Q = W H into (36a) and let Z = PP H . The result is V = APR for 
some unitary R. Therefore, the linear system of (42) has a solution Uo = PR- It is easy to see that /(Uo) = 
I(P) = /(Z), thus (36d) results in |/(Uo)| < <r(A)/2. Applying Theorem 2 to each of the columns of Uo 
provides the uniqueness of Uo- It is trivial that Z = UoU^. ■ 

Using the same arguments as in the proof it is easy to conclude that I(Z) = /(Uo), so that St can be found 
directly from the solution matrix Uo- In particular, we develop the Continuous to Finite (CTF) block which 
determines the diversity set St of a given frequency interval T. Fig. 2 presents the CTF block that contains the 
flow of transforming the continuous linear system of (16) on the interval T into the finite dimensional problem 
of (42) and then to the recovery of St- The role of Propositions 2 and 3 is also illustrated. The CTF block is the 
heart of the SBR scheme which we discuss next. 

In the CS literature, the linear system of (42) is referred to as an MMV system. Theoretical results regarding 
the sparsest solution matrix of an MMV system are given in [15]. Finding the solution matrix Uo is known to 
be NP-hard [12]. Several sub-optimal efficient algorithms for finding Uo are given in [16]. Some of them can 
indicate a success recovery of Uq. We explain which class of algorithms has this property in Section VI-A. 
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ylf) 

V/6 T 

T . 



Q = Jy(f)y H (f)df 

T 



Q = W" 



MMV problem 



AU 



Solve MMV 

Brute-force or using 
CS algorithms 



S T = J(U ) 



St 



Set of constraints (36) 
Unique solution 



Unique solution 
Proposition 3 



Theoretical 
results 



Proposition 2 



Fig. 2. Continuous to finite block (CTF). This block determines the diversity set St of a given interval T. 



VI. SBR ALGORITHMS 

The theoretical results developed in the previous section are now used in order to construct the diversity set S 
which enables the recovery of x(t) via (29)-(30). 

We begin by defining a class A of signals. The SBR4 algorithm is then presented and is proved to guarantee 
perfect reconstruction for signals in A. We then show that in order to ensure that M. C A the sampling rate must 
be at least 4NB, which is twice the minimal rate stated in Theorem 1. To improve on this result, we define a 
class B of signals, and introduce a conceptual method to perfectly reconstruct this class. The SBR2 algorithm is 
developed so that it ensures exact recovery for a subset of B. We then prove that M. is contained in this subset 
even for sampling at the minimal rate. However, the computational complexity of SBR2 is higher than that of 
SBR4. Since universal patterns lead to the largest sets A and B, we assume throughout this section that universal 
patterns are used, which results in <r(A) = p. 

A. The SBR4 algorithm 

Define the class Ak of signals 

A K = {suppX(/) C T and \S\ < K}, (43) 

with S given by (25). Let T = Fq, and observe that a multi-coset system with p > 2K ensures that all the 
conditions of Proposition 2 are valid for every x(t) G Ak- Thus, applying the CTF block on T = JT results in 
a unique sparsest solution Uo, with S = I(Uo). The reconstruction of the signal is then carried out by (29)-(30). 
We note that (27) is valid as it represents the class A p that contains Ak for p > 2K. 

Algorithm 1, named SBR4, follows the steps of the CTF block in Fig. 2 to recover the diversity set S from 
y(/), for any x(t) G Ak- The algorithm also outputs an indication flag which we discuss later on. 

The SBR4 algorithm guarantees perfect reconstruction of signals in Ak from samples at twice the Landau rate, 
which is also the lower bound stated in Theorem 1. To see this, observe that (25) implies that every x(t) G Ak 
must satisfy 

A(suppX(/)) < A (44) 
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Algorithm 1 SBR4 
Input: y(/), Assume: a (A) = p 
Output: the set S, flag 
l: Set T = 

2: Compute the matrix Q by (33) 

3: Decompose Q = VV^ according to Proposition 3 

4: Solve the MMV system V = AU for the sparsest solution Uo 

5: S = /(Uo) 

6: flag={|S|<2} 

7: return 5, flag 



Although Ak is not a subspace, we use (44) to say that the Landau rate for Ak is K/LT as it contains subspaces 
whose widest support is K/LT. As we proved, p > 2K ensures perfect reconstruction for Ak- Substituting the 
smallest possible value p = 2K into (19) results in an average sampling rate of 2K/LT. 

It is easy to see that flag is equal to 1 for every signal in Ak- However, when a sub-optimal algorithm is used 
to solve the MMV in step 4 we cannot guarantee a correct solution Uo- Thus, flag=0 indicates that the particular 
MMV method we used failed, and we may try a different MMV approach. 

Existing algorithms for MMV systems can be classified into two groups. The first group contains algorithms 
that seek the sparsest solution matrix Uo, e.g. Basis Pursuit [17] or Matching Pursuit [18] with a termination 
criterion based on the residual. The other contains methods that approximate a sparse solution according to user 
specification, e.g. Matching Pursuit with a predetermined number of iterations. Using a technique from the latter 
group neutralizes the indication flag as the approximation is always sparse. Therefore, this set of algorithms should 
be avoided if an indication is desired. 

An important advantage of algorithm SBR4 is that the matrix Q can be computed in the time domain from the 
known sequences x Ci [n], 1 < i < p. The computation involves a set of digital filters that do not depend on the 
signal and thus can be designed in advance. The exact details are given in Appendix C. 

The drawback of the set Ak, is that typically we do not know the value of K. Moreover, even if K is known, 
then usually we do not know in advance whether x(t) 6 Ak as Ak does not characterize the signals according 
to the number of bands and their widths. Therefore, we would like to determine conditions that ensure M. C Ak- 
Proposition 1 shows that for x(t) G M. the set S satisfies \S\ < 2N if L < 1/BT. Thus, under this condition on 
L we have M. C A2N, which in turn implies p = AN as a minimal value for p. Consequently, SBR4 guarantees 
perfect reconstruction for M. under the restrictions L < 1/ BT and p > AN. However, the Landau rate for M. is 
NB, while p = AN implies a minimal sampling rate of ANB. Indeed, substituting p = AN and L < 1/BT into 
(19) we have 

p AN 

— > — =- = ANB. (45) 
LT ~ T-^ K ' 

In contrast, it follows from Theorem 3 that p > 2N is sufficient for uniqueness of the solution. The reason for 
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the factor of two in the sampling rate is that x(/) is iV-sparse for each specific /; however, when combining the 
frequencies, the maximal size of S is 2N. The SBR2 algorithm, developed in the next section, capitalizes on this 
difference to regain the factor of two in the sampling rate, and thus achieves the minimal rate, at the expense of 
a more complicated reconstruction method. 

B. The SBR2 algorithm 

We now would like to reduce the sampling rate required for signals of M. to its minimum, i.e. twice the Landau 
rate. To this end, we introduce a set Bk for which SBR2 guarantees perfect reconstruction, and then prove that 
MCB N if p> 2N. 

Consider a partition of JF into M consecutive intervals defined by 

= di < d 2 < ■ ■ ■ < d~M+i = j7p- 
For a given partition set D = {di} we define the set of signals 

B K ,D = {snppX(f) C T and < K, 1 < i < M}. 

Clearly, if p > 2K then we can perfectly reconstruct every x(t) € B K ^ by applying the CTF block to each of 
the intervals [di,di + \\. We now define the set Bk as 

Bk = \JB k>3 , (46) 

D 

which is the union of B K q over all choices of partition sets D and integers M. Note that neither Bk nor B K q 
is a subspace. If we are able to find a partition D such that x(t) € B K q, then x(t) can be perfectly reconstructed 
using p > 2K. Since the Landau rate for Bk is K/LT, this approach requires the minimal sampling rate 2 . 

The following proposition shows that if the parameters are chosen properly, then M. C Bn- Thus, p > 2N and 
a method to find D of x(t) is sufficient for perfect reconstruction of x(t) € M.. 

Proposition 4: If L,p,C are selected according to Theorem 3 then M. C Bn- 

Proof: In the proof of Theorem 3 we showed that under the conditions of the theorem, x(/) is iV-sparse for 
every / € Tq. The proof of the proposition then follows from the following lemma [8]: 

Lemma 2: If x(t) is a multi-band signal with N bands sampled by a multi-coset system then there exists a 
partition set D = {di} with M = 2N + 1 intervals such that I(x(/)) is a constant set over the interval [d~i,d~i + i] 
for 1 < % < M. 

Lemma 2 implies that \Sa. j. i j < N for every 1 < i < M = 2N + 1 which means that x(t) £ B N q. ■ 
So far we showed that M. C Bn, however to recover x(t) we need a method to find D in practice; Lemma 2 



2 under the convention discussed for Ak- 
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only ensures its existence. Given the data y(/), our strategy is aimed at finding any partition set D such that 

\D\-1 

S= |J S [diA+l] (47) 

i=0 

is equal to S, and such that l^r^di+ill — ^ ^ or ever y 1 < ^ < M. As long as (27) holds, once we find S 
the solution is exactly recovered via (29)-(30). To find S, we apply the CTF block on each interval [di,di + i]. If 
p > 2K, then the conditions of Proposition 2 are valid, a unique solution is guaranteed for each interval. Since for 
p = 2K (27) is valid for A 2 k, our method guarantees perfect reconstruction of signals in Bk H A 2 k- As always, 
using a universal pattern makes the set of signals Bk H A 2 k the largest. Since the Landau rate for Bk H A 2 k is 
K/LT this approach allows for the minimal sampling rate when p = 2K. 

In order to find D we suggest a bi-section process on Tq. We initialize T = Tq and seek St- If Sj- does 
not satisfy some condition explained below, then we halve T into T\ and T 2 and determine Sq\ and SV 2 . The 
bi-section process is repeated several times until the conditions are met, or until it reaches an interval width of no 
more than e. The set S is then determined according to (47). 

We now describe the conditions for which a given T C JF is halved. The matrix Z of (34) satisfies the 
constraints (36a)-(36b). Since x(t) € A2K and p > 2K (36c) is also valid. However, the last constraint (36d) of 
Proposition 2 is not guaranteed as it requires a stronger condition \St\ < K = p/2. Note that this condition is 
satisfied immediately if D = D since x(t) E Bk- We suggest to approximate the value Sj- = \I(Z ) \ by rank(Q), 
and solve the MMV system for the sparsest solution only if rank(Q) < p/2. This approximation is motivated by 
the fact that for any Z y it is true that rank(Z) < |/(Z)|. From Proposition 2 we have that rank(Z ) = rank(Q) 
which results in 

rank(Q) < |/(Z)|. (48) 

However, only special multi-band signals result in strict inequality in (48). Therefore, an interval T that produces 
rank(Q) > p/2 is halved. Otherwise, we apply the CTF block for this T assuming that (48) holds with equality. 
As in SBR4 the flag indicates a correct solution for x(t) G Bk fl A 2 k- Therefore, if the flag is we halve T. 
These reconstruction steps are detailed in Algorithm 2, named SBR2. 

It is important to note that SVR2 is sub-optimal, since the final output of the algorithm S may not be equal to 
S even for x(t) 6 Bk H A 2 k- One reason this can happen is if strict inequality holds in (48) for some interval 
T. In this scenario step 7 is executed even though Zo does not satisfy (36d). For example, a signal x{t) with two 
equal width bands [01, a\ + W] and [a 2 , a 2 + W\ such that 



«1 




d2 


ILT- 







and 7 + W € Tq. If x(t) also satisfies 

X(f- ai ) = X(f-a 2 ), V/E[(W], (50) 
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Algorithm 2 SBR2 

Input: T, Initialize: T = Tq, Assume: cr(A) = p 
Output: a set S 

l: if A(T) < e then 

2: return 5 = {} 

3: end if 

4: Compute the matrix Q by (33) 

5: if rank(Q) < \ then 

6: Decompose Q = W H 

l: Solve MMV system V = AU 

8: S = 7(U ) 

9: else 

10: 5 = {} 

11: end if 

12: if (rank(Q) > f ) or (\S\ > §) then 

13: split T into two equal width intervals T\, T2 

14: £W = SBR2(Ti) 

15: §W = SBR2(T 2 ) 

16: S = U 

17: end if 
18: return S 



then it can be verified that |J(Zo)| = 2 while rank(Zo) = rank(Q) = 1 on the interval T = [7,7 + W]. This is of 
course a rare special case. Another reason is a signal for which the algorithm reached the termination step 1 for 
some small enough interval. This scenario can happen if two or more points of D reside in an interval width of e. 
As an empty set S is returned for this interval, the final output may be missing some of the elements of S. Clearly, 
the value of e influences the amount of cases of this type. We note that since we do not rely on D = D the missing 
values are typically recovered from other intervals. Thus, both of these sources of error are very uncommon. 

The most common case in which SBR2 can fail is due to the use of sub-optimal algorithms to find Uo; this 
issue also occurs in SBR4. As explained before, we assume that fiag=0 means an incorrect solution and halves the 
interval T. An interesting behavior of MMV methods is that even if Uo cannot be found for T, the algorithm may 
still find a sparse solution for each of its subsections. Thus, the indication flag is also a way to partially overcome 
the practical limitations of MMV techniques. Note that the indication property is crucial for SBR2 as it helps to 
refine the partition D and reduce the sub-optimality resulting from the MMV algorithm. 

We point out that Proposition 4 shows that M. C Bn- We also have that M. C A2N from Proposition 1, 
which motivates our approach. The SBR2 algorithm itself does not impose any additional limitations on L,p,C 
other than those of Theorem 3 required to ensure the uniqueness of the solution. Therefore, theoretically, perfect 
reconstruction for A4 is guaranteed if the samples are acquired at the minimal rate, with the exception of the 
special cases discussed before. 

The complexity of SBR2 is dictated by the number of iterations of the bi-section process, which is also affected 



TABLE II 

Spectrum-blind reconstruction methods for multi-band signals 





WKS theorem 


SBR4 


SBR2 


Sampling method 


Uniform 


Multi-coset 


Multi-coset 


Fully-blind 


Yes 


Yes 


Yes 


# Uniform sequences 


1 


P 


P 


Minimal sampling rate 


Nyquist 


2 x Landau 


2 x Landau 


Achieves lower bound of Theorem 1 


No 


Yes 


Yes 


Reconstruction method 


Ideal low pass 


SBR4 


SBR2 


Time complexity 


constant 


1 MMV system 


bi-section + finite # of MMV 


Applicability 


suppX(/) C T 


x(t) e Ak 


x(t) e B K n A 2K 3 


Indication 


No 


for x(t) <G Ak only 


No 



by the behavior of the MMV algorithm that is used. Numerical experiments in Section VII show that empirically 
SBR2 converges sufficiently fast for practical usage. 

Finally, we emphasize that SBR2 does not provide an indication on the success recovery of x(t) even for 
x(t) € M since there is no way to know in advance if x(t) is a signal of the special type that SBR2 cannot 
recover. 

C. Comparison between SBR4 and SBR2 

Table II compares the properties of SBR4 and SBR2. We added the WKS theorem as it also offers spectrum- 
blind reconstruction. Both SBR4 and SBR2 algorithms recover the set S according to the paradigm stated in 
Section V-B. Observe that an indication property is available only for SBR4 and only if the signals are known to 
lie in Ak- Although both SBR4 and SBR2 can operate at the minimal sampling rate, SBR2 guarantees perfect 
reconstruction for a wider set of signals as Ak is a true subset of Bk H A%k- 

Considering signals from M. we have to restrict the parameter selection. The specific behavior of SBR4 and 
SBR2 for this scenario is compared in Table III. In particular, SBR4 requires twice the minimal rate. 

In the tables, perfect reconstruction refers to reconstruction with a brute-force MMV method that finds the 
correct solution. In practice, sub-optimal MMV algorithms may result in failure of recovery even when the other 
requirements are met. The indication flag is intended to discover these cases. 

The entire reconstruction scheme is presented in Fig. 3. The scheme together with the tables allow for a wise 
decision on the particular implementation of the system. Clearly, for > 0.5 it should be preferred to sample at 
the Nyquist rate and to reconstruct with an ideal low pass filter. For < 0.5 we have to choose between SBR4 
and SBR2 according to our prior on the signal. Typically, it is natural to assume x(t) G M for some values of 



3 except for special signals discussed in Section VI-B. 



TABLE III 

Comparison of SBR4 and SBR2 for signals in M 
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SBR4 


SBR2 


# Uniform sequences 
Minimal rate 
Lower bound of Th. 1 
Parameter selection 
Perfect reconstruction 
Indication 


p > AN 
4 x Landau 
No 

Theorem 3, p > AN 
Yes 
Yes 


p > 2N 
2 x Landau 

Yes 
Theorem 3 
Yes 3 
No 




Fig. 3. Spectrum-blind reconstruction scheme. 



N and B and derive the required parameter selection according to Table III. It is obvious that if p > AN is used 
then SBR4 should be preferred since it is less complicated than SBR2. 

The trade-off presented here between complexity and sampling rate also exists in the known-spectrum reconstruc- 
tion of [8]. Sampling at the minimal rate of Landau requires a reconstruction that consists of piecewise constant 
filters. The number of pieces and the reconstruction complexity grow with L. This complexity can be prevented 
by doubling the value of p which also doubles the average sampling rate according to (19). Then, (29)-(30) are 
used to reconstruct the signal by only one inversion of a known matrix [6]. 

VII. Numerical experiments 

We now provide several experiments demonstrating the reconstruction using algorithms SBR4 and SBR2 for 
signals from M.. We also provide an example in which the signals do not lie in the class M. but in the larger set 
implied by A K for SBR4 and by B K n A 2 k for SBR2. 
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A. Setup 

The setup described hereafter is used as a basis for all the experiments. 

Consider an example of the class M with T = [0,20 GHz],iV = 4 and B = 100 MHz. In order to test the 
algorithms 1000 test cases from this class were generated randomly according to the following steps: 

1) draw {a{}fL 1 uniformly at random from [0, 20GHz — B}. 

2) set bi = cii + B for 1 < i < N, and ensure that the bands do not overlap. 

3) Generate X(f) by 



N 

«(/)(£«(/) /e UK A] 

i=l 



, 0, otherwise. 
For every / the values of Sn(f) and Sj(f) are drawn independently from a normal distribution with zero 
mean and unit variance. The function a(f) is constant in each band, and is chosen such that the band energy 
is equal to where e.- t is selected uniformly from [1,5]. 

The Landau rate for each of the signals is NB = 400 MHz, and thus the minimal rate requirement for blind 

reconstruction is 800 MHz due to Theorem 1. 

Several multi-coset systems are considered with the following parameters. The value L is common in all the 

systems. The value of p is varied from p = N = A to p = 8N = 32 representing 29 different systems. A universal 

pattern C is constructed by choosing prime L, since according to [10] this ensures that every sampling pattern is 

universal. 

An experiment is conducted by sampling the signals using each of the multi-coset systems. Each of these 
combinations is used as an input to both SBR4 and SBR2 algorithms. We selected the Multi-Orthogonal Matching 
Pursuit (M-OMP) method [16] to solve the MMV systems for the sparsest solution. The empirical success rate of 
each algorithm is calculated as the ratio of simulations in which the recovered set S is correct. 

B. Sampling rate and practical limitations 

We begin by selecting the largest possible value of prime L satisfying (23): 

L = 199 < = 200. (51) 

Thus, the minimal rate requirement holds only for p > 2N. Specifically, for p = 2N the sampling rate is 
p/LT = 804 MHz. Observe that a non-prime L = 200 would give the minimal rate exactly. This setting is 
discussed later on. 

Fig. 4 depicts the empirical success rate with L = 199, N = 4 as a function of p. It is evident that for p < 2N the 
set S could not be recovered by neither of the algorithms since the sampling rate is below the bound of Theorem 1 . 
As expected, SBR2 outperforms SBR4 as it achieves the same empirical success rate for a lower average sampling 
rate. It is also seen that for p = 4N the sampling rate is slightly more than four times the Landau rate. Indeed, 



24 




Fig. 4. Performance of SBR algorithms with L — 199. 
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Fig. 5. Performance of SBR algorithms with L — 23. 



algorithm SBR4 maintains a high recovery rate for this value of p. The usage of SBR2 with M-OMP maintains a 
high recovery rate for p/N = 2.6, which is more than the minimal rate. Other MMV algorithms may be used to 
improve this result, however we used only M-OMP as it is simple and fast. 

We next consider a scenario with L = 23, which clearly satisfies (23). Here, for p = N = A we have a sampling 
rate of 3.4 GHz which is much higher than the minimal requirement. This selection of L represents a practical 
desire to satisfy the minimal rate requirement with a reduced value of p, since realizing the multi-coset sampling 
requires p analog-to-digital devices. Fig. 5 presents the empirical recovery rate in this case. Note that Table III 
shows that in order to guarantee perfect reconstruction for M. we need p > AN for SBR4, and p > 2N for 
SBR2. However, these conditions are only sufficient. Indeed, it is evident from Fig. 5 that both algorithms reach 
a satisfactory recovery rate for lower values of p . 

In Table IV, we tabulate the average run time of one case out of the 1000 tested. Our experiments were conducted 
on an ordinary PC desktop with an Intel CPU running at 2.4GHz and 512MB memory RAM. We used Matlab 
version 7 to encode and execute the algorithms. Note that for L = 199, p = 2N we encountered a significant 
increase in SBR2 runtime. The reason is that the average sampling rate is very close to the minimal possible, 
thus the recursion depth of the algorithm grows as it is harder to find a suitable partition set D. For p = AN the 



TABLE IV 

Average run time of SBR4 and SBR2 with MOMP (msec) 
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L = 


199 


L = 


: 23 


SBR4 


SBR2 


SBR4 


SBR2 


p = 


N 


7 


608 


4.2 


51.4 


p = 


2N 


16.1 


1034 


5.7 


6.4 


p = 


AN 


21.4 


24.8 


6.7 


6.7 



runtime dramatically improves, however in this case SBR4 may be preferred due to the advantages that appear in 
Table III. It can be seen that for L = 23 the average runtime is low for both algorithms. This scenario represents 
a case that the value of |5| is very low compared to 2N, and thus it is easier to find a partition set D. Moreover, 
M-OMP becomes faster as the solution is sparser. 

C. Applicability 

The previous experiments demonstrated the applicability of SBR4 and SBR2 to signals that lie in M. We now 
explore the case in which x(t) £ Ai. 

In this experiment we used the basic setup with L = 199 but the signals are constructed in a different way. Each 
one of the 1000 signals is constructed by X(f) = a(f) (»Sr(/) + jSi(f)) , V/ G Fq. The function a(f) depends 
on the algorithm and it makes sure that x(t) G Ak for the test cases of SBR4. Similarly, a(f) is used to form 
signals x(t) G Bk H A%k for SBR2. The construction of these signals depends on L because of the definitions of 
Ak an d Bk- We selected K = 8 which results in a Landau rate of K/LT = 804 MHz in either construction. In 
addition, we made sure that the signals do not lie in M.. 

Fig. 6 shows the empirical recovery rate of SBR4 and SBR2 in this scenario. The value p = AN = 16 serves as 
a threshold for satisfactory recovery, as the sampling rate for this value of p is p/LT = 1608 MHz, which is twice 
the Landau rate. It can also be seen that SBR4 performs better than SBR2 as it does not involve a sub-optimal stage 
of recovering the partition set D. Both algorithms suffer from the sub-optimality techniques for MMV systems. 

Note that the signals here are synthesized so that they lie in the relevant sets. However, for a generic signal 
x(t) £ M. there is no way to know in advance whether it lies in one of these sets. Moreover, there is no way to 
infer it from the samples, y(/). In addition, even if SBR4 is used for this signal and it returns flag=l, there is no 
meaning for this indication since the uniqueness of the solution is guaranteed only for x(t) G Ak which cannot 
be ensured for a generic multi-band signal. 

D. Random sampling patterns 

Theorem 3 requires a universal sampling pattern, which means finding a pattern resulting in <r(A) = p. However, 
computing the value of a (A) requires a combinatorial process for non-prime L. The "bunched" pattern C = 
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Fig. 6. Performance for signals x(t) not in M. 




Fig. 7. Performance of SBR algorithms with L — 200. 



{0, 1, ...,p — 1} given in [8] is proved to be universal but the matrix A is not well conditioned for this choice 
[7]. Alternatively, it follows from the work of Candes et. al. [19] that random sampling patterns are most likely 
to produce a high value for a (A) if L,p are large enough. Therefore, for practical usage the sampling pattern 
can be selected randomly even for non-prime L. Fig. 7 presents an experiment with L = 200. We point out that 
the random selection process is carried out only once, and the same sampling patterns are used for all the tested 
signals. Comparing Figs. 4 and 7 it is seen that the results are very similar although the exact value of a (A) is 
unknown. 

This experiment was also performed when for every N < p < 8N the patterns are selected as C = {ck\ck = 
2k, < k < p — 1}, which is proved in [10] to render a (A) = 1. In this case both SBR4 and SBR2 could not 
recover any of the 1000 test cases. Thus, the universality of the pattern is crucial to the success of our method. 

VIII. Conclusions 

In this paper we suggested a method to reconstruct a multi-band signal from its samples when the band 
locations are unknown. Our development enables a fully spectrum-blind system where both the sampling and 
the reconstruction stages do not require this knowledge. 
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Our main contribution is in proving that the reconstruction problem can be formulated as a finite dimensional 
problem within the framework of compressed sensing. This result is accomplished without any discretization. 
Conditions for uniqueness of the solution and algorithms to find it were developed based on known theoretical 
results and algorithms from the CS literature. 

In addition, we proved a lower bound on the sampling rate that improves on the Landau rate for the case of 
spectrum-blind reconstruction. One of the algorithms we proposed indeed approaches this minimal rate for a wide 
class of multi-band signals characterized by the number of bands and their widths. 

Numerical experiments demonstrated the trade off between the average sampling rate and the empirical success 
rate of the reconstruction. 



In order to treat real-valued signal the following definitions replace the ones given in the paper. The class M 
is changed to contain all real-valued multi-band signals bandlimited to T = [-1/2T, 1/2T] with no more than N 
bands on both sides of the spectrum, where each the band width is upper bounded by B as before. Note that N is 
even as the Fourier transform is conjugate symmetric for real- valued signals. The Nyquist rate remains 1/T and 
the Landau rate is NB. 

Repeating the calculations of [8] that lead to (15) it can be seen that several modifications are required as 
now explained. To form x(/), the interval T is still divided into L equal intervals. However, a slightly different 
treatment is given for odd and even values of L, because of the negative side of the spectrum. Define the set of 
L consecutive integeres 



Appendix A 



Real-valued signals 




odd L 



even L. 



and redefine the interval JT 



1 



1 



odd L 



2LT 2LT 



= < 




even L. 



The vector x(/) is now defined as 



Xi(f)=X(f + Ki/LT), V0<i<L-l, 



2S 



The dimensions of A remain p x L with ik entry 

A ' k = Jjf exp ( 3 T CiKt ) ' (52) 

l<i<p, < k < L- 1. 

The definition of y(/) remains the same with respect to JF defined here. The results of the paper are thus extended 
to real-valued multi-band signals since (16) is now valid with respect to these definitions of x(/), A, and JF . 

Note that, we could have, conceptually, constructed a complex-valued multi-band signal by taking only the 
positive frequencies of the real- valued signal. The Landau rate of this complex version is NB/2. Nevertheless, 
the information rate is the same as each sample of a complex-valued signal is represented by two real numbers. 

Appendix B 
Proof of Lemma 1 

Let r = rank(P). Reorder the columns of P so that the first r columns are linearly independent. This operation 
does not change the rank of P nor the rank of AP. Define 

P = [P (1) P (2) ], (53) 

where P^ 1 ) contains the first r columns of P and the rest are contained in P^ 2 ). Therefore, 

r > rank(AP) = rank(A[p( 1 ) p( 2 )]) > rank(AP( 1) ). 

The inequalities result from the properties of the rank of concatenation and of multiplication of matrices. So it is 
sufficient to prove that AP^ has full column rank. 

Let a be a vector of coefficients so that AP^a = 0. It remains to prove that this implies a = 0. Denote 
k = |/(P)|. Since I(P^) C I(P) = k the vector P^a is A>sparse. However, a(A) > k and its null space 
cannot contain a A;-sparse vector unless it is the zero vector. Since P^ 1 ) contains linearly independent columns this 
implies a = 0. 

Appendix C 
Computation of the matrix Q 

The SBR4 algorithm computes the matrix Q in the frequency domain. A method to compute this matrix directly 
from the samples in the time domain is now presented. 
Consider the ikth element of Q from (33): 

Q ik = f LT Yi(f)yUf)df. (54) 
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Since yi(f) is the DTFT of x Ci [n] we can write Qi k as, 



Qik = / [yz Xc - [ n *] ex p (-j 27r f n i T ) • ( 55 ) 

^ Wz / 

I ^2 x Ck [n k ]exp(-j2irfn k T) J df 

1 

^ ^2 x c, [ni]x* Ck [n k ] / exp (j2vr/(n fc - n^T) (if. 



Note that from (14) the sequence x Ci [nj] is padded by L — 1 zeros between the non-zero samples. Define the 
sequence without these zeros as 

x Ci [m] = x(mLT + qT), m £ Z, 1 < i < p. (56) 

Then, (55) can be written as 

Qik=^2 ^2 x Ct [mi\x* Ck [m k ]g ik [mi - m k ] (57) 



where 



and 



LT 



aik[m}= exp(j2irf(mL + (c k - Ci))T)df, (58) 



n 



\x ru * 



9ik) M = ^2 ^ 9ik I m ~ n ] • ( 59 ) 



If i = k then Cj = and 

= #[m] = --^ exp(j7rm) sinc(m), (60) 

with sinc(x) = sin(-7rx)/(7ra;). 
If i ^ k, 

I l ^p (j2"( Ck -Ci))-l 

9ik[m\ = . (61) 

j2vr(mL + (c k - Ci))T 

The set of digital filter g ik can be designed immediately after setting the parameter L,p,C as these filters do 
not depend on the signal. 
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